Monte-Carlo rejection as a tool for measuring the energy landscape scaling of simple 

fluids 
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A simple modification of the Monte-Carlo algortihm is proposed to explore the topography and 
the scaling of the energy landscape. We apply this idea to a simple hard-core fluid. The results for 
different packing fractions show a power law scaling of the landscape boundary, with a characteristic 
scale that separates the values of the scaling exponents. Finally, it is shown how the topology 
determines the freezing point of the system due to the increasing importance and complexity of the 
boundary. 
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I. INTRODUCTION 



A liquid cooled to temperatures near its freezing point 
can be conduced to a glassy state or to a crystal accord- 
ing to the speed of cooling |lj. When the speed is high 
enough, the supercooled liquid undergoes a glass tran- 
sition to a state that is disordered, while a phase tran- 
sition of the type fluid-crystal is obtained if the system 
is kept in thermodynamical equilibrium at all steps of 
the cooling path 2] . The understanding of the many dif- 
ferent aspects of glass transition still remains as one of 
the most important problems in physics as for ex- 
ample, the explanation of the non-exponential relaxation 
of fluctuations £| or why not all materials form glasses 
5j. Another very interesting property of glasses, related 
with the glass forming tendency jo] , is the behavior of the 
viscosity, which is usually referred as fragility. Different 
approaches have been used to understand glass transi- 
tion : models like the Gibbs-DiMarzio [jj, theories like 
the mode coupling, stochastic agglomeration [HIHEl or 
the use of computer simulations 0. Another useful ap- 
proach is the rigidity theory of Phillips and Thorpe 
[L2j . which relates the ratio between available degrees of 
freedom and the number of constraints [13j. In previous 
works, we showed that even for simple systems like hard- 
disks 0| and colloids [l^ , it seems that rigidity plays an 
important role even for the case of a simple phase tran- 
sition, since it is clear that in order to form a solid, the 
system must develop certain rigidity. Some works on the 
relaxation properties of colloids, seem to confirm these 
ideas [pq . 

Parallel to all of these approaches, there is another for- 
malism that has been very useful to visualize and under- 
stand what happens during a glass or usual phase tran- 
sition . This formalism is the energy landscape approach 
|17| [Tsj , which many years ago was very successful in the 
field of spin glasses Jlj| ■ The main idea behind this ap- 
proach, is that the landscape is a surface generated by the 
potential energy of the system as a function of the molec- 
ular positions For a system with N molecules, the 
landscape is determined by the potential energy function, 



<I>(ri, ...,rjv), where r-j comprise all relevant coordinates, 
like position and orientation. Since the kinetic energy 
(K) is a positive defined quantity, the system evolves in 
such a way that K = E — ^(ri, rjy) > 0, where E 
is the total energy. The topography of the landscape 
energy surface determines the possible evolution of the 
system, and the contact with thermodynamics is made 
by using statistical mechanics 0- The great advantage 
of the landscape formalism is that it can be used even 
without thermodynamical equilibrium. In such a case, 
ergodicity is broken and the entropy is not a maximum 
anymore, as postulated in the usual equilibrium statis- 
tical mechanics. The main question to be answered for 
this case, is to figure out the fraction of allowed volume 
in phase space that is visited by the system. 

The usual picture of a phase or glass transition in such 
a language, is that at high temperatures the system does 
not feel the topology of 3>(ri, rjy) because the kinetic 
energy contribution dominates. As the temperature is 
lowered, the system is unable to surmount the highest 
energy barriers and therefore is forced to sample deep 
minima. For a slow cooling is slow, the system has time 
to find a path to the most stable state, an ordered crys- 
tal. It will be trapped in a metastable state, the glass, 
[l7| fl8| if the cooling speed is high enough. Many works 
that relates the statistics of an energy landscape with 
the thermodynami cal p roperties of the system have been 
made [2(| |2l| |23|. and even some phenomenological 
relations between rigidity and the energy landscape have 
been obtained (24[. For a Lennard- Jones fluid, it has 
been determined that the network of minima is a static 
scale- free network pB|. 

However, up to now there are some important ques- 
tions that remain, for example, what is the topography 
of the landscape for a fragile liquid or how to obtain the 
viscosity jlSj|. diffusion |2q| and rigidity from the land- 
scape [24J. Also, although is widely believed that the 
landscape is fundamental to understand many features 
of a liquid, still is not completely clear how to use the 
landscape topology to predict a phase transition p^ . 
Another interesting question is: what is the nature of 
the texture of the landscape? In other words, do the to- 
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pography is a fractal? What's the fractal dimension of 
the landscape? Although some of these questions seem 
to have an academic interest only, is clear that the relax- 
ation properties of a very complex fractal landscape are 
different from a smooth landscape , where the system 
can easily explore the phase space. In this article, we will 
explore some of these questions by looking at the scaling 
of the landscape. 

The layout of this work is the following: in section II 
we discuss how to use a modified Monte-Carlo method 
to study the scaling. In section III we apply the method 
to a simple fluid of hard disks. Finally, in section IV we 
give the conclusions. 

II. MONTE-CARLO REJECTION AND 
SCALING OF THE LANDSCAPE 

In this section we will develop a method to relate the 
Monte-Carlo rejection ratio and the scaling of the land- 
scape. Before going into the details, it would be useful 
to explain some others approaches to obtain information 
about the landscape topology. To simplify ideas, in this 
article we will use as a model system, N hard disks or 
spheres of diameter a in a given area (S) or volume ( V) . 
In such hard-core particle system, the energy landscape 
is formed by walls of infinite height that divide the al- 
lowed and forbidden regions of the configurational phase 
space. If Yi is the position of a disk or sphere i, the al- 
lowed region of the landscape is the set of points where, 

ll^-r.H >a, (1) 

for all possible pairs i and j. The number of such pairs 
(Rnh) is the number of combinations of N objects taken 
in pairs: = N(N — l)/2. Each of these equations 
is a non-holonomic restriction to the system. A state P 
in phase space is in the boundary of the landscape, if at 
least one of the inequalities is transformed into an 
holonomic restriction, 

||rj - r m || = a, (2) 

for a pairs of disks that we denote by I and m. For each 
equation of this type that is satisfied, two disks are in 
contact. For a given packing fraction (</>), the number of 
such equations (Rh) is just, 

«. - (3) 

where (Z (</))) is the average coordination per particle in 
a given packing. We remark that this equation allows 
a straight forward manner to connect the energy land- 
scape formalism with the rigidity theory, where the most 
important parameter is {Z ((/>)). This approach also 
provides a way to construct inherent structures and the 
boundary of the landscape just by considering a non- 
linear optimization problem. To get a packing, we can 



define an objective function as, 

JY 

F(P) = X>*II. (4) 

3=1 

with Rh constraints ||r; — r m || = er. This objective func- 
tion is defined in such a way that the particles are packed 
in a tight way with respect to the origin. Surprisingly, 
this criteria is similar to the center of mass minimization 
that has been observed very recently in experiments with 
colloids 29] 30] . In principle, this problem can be solved 
using nonlinear programming [3l| , and there is some sim- 
ilar work done into this line of research |32| . Here we 
will not follow this path. Instead, we investigate how the 
landscape boundary looks at different scales in the config- 
urational part of the phase space. The most simple way 
to do this, consists in applying a box-counting algorithm 
[33] ] once the boundary points are determined. In this box 
counting algorithm, a grid made from cubes of linear size 
5 is applied to the configurational phase space. Then, 
the number of boxes that contains a boundary state are 
counted. The counting is repeated at different lengths 
5. This ideal situation has the problem that we need 
to generate all the boundary states, and due to the size 
of the phase space, this task is almost impossible to do. 
A simpler approach is to take advantage of the Monte- 
Carlo importance sampling to obtain information about 
the boundary. 

In a usual Metropolis Monte-Carlo j^J 35], when a 
system is in a microstate Po, a movement to a new mi- 
crostatc P is allowed if the difference in energies AE = 
E(P) — E(Pq) is negative (where E(P) is the energy of 
the state P). If AE > 0, a random number is compared 
with AE. In the case of hard-core systems, a rejection 
occurs when a new proposed configuration is not allowed 
by the restrictions |35|. i.e., when there is an overlap be- 
tween disks or spheres. When the new proposed point 
P is rejected, one can be sure that the boundary of the 
landscape has been crossed, i.e., the boundary is between 
states Po and P. Thus, the information about the bound- 
ary can be extracted from the acceptance ratio of the 
Monte-Carlo, although two modifications are needed. In 
the usual Monte-Carlo, the trial movement distance be- 
tween two states is a continuous random variable |36| . 
This feature is not convenient because it does not pro- 
vide an approximate location for the boundary. A second 
fact to take into account, is that the probability of hitting 
the point P not only depends on the size of the boundary, 
but also in the transient probability of the process j37| 
Po — > P. To solve these problems, let us introduce a reg- 
ular grid in the configurational part of the phase space. 
If the simulation is performed in a box of linear length 
L, there are M = (L/5) DN points in the grid, where D 
is the dimensionality of the system. In such a grid, a 
random walk in phase space is performed by choosing at 
random a disk and changing one components at ran- 
dom to ±8. If P is written in generalized coordinates qj, 
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the trial movement is written as, 

(QU ?2,«3, Qdn) -> (qi, 92, 93, q r ± <5, tow), (5) 

for the r coordinate chosen at random with an uniform 
distribution. 

The introduction of a random walk is convenient be- 
cause: 1) the step size can be varied to look at the scaling 
and 2) when a movement is rejected, the state can be con- 
sidered as a boundary point, since it is connected to the 
interior of the forbidden part of the phase space. In spite 
of this, the introduction of a random walk has the incon- 
venient of not being able to sample the phase space with 
equal probabilities, i.e., it is not completely ergodic. In 
fact, even for the usual Monte-Carlo there is some lack of 
ergodicity due to the existence of an finite-size underly- 
ing grid. The only difference between the simple random 
walk and the usual Monte-Carlo with jumps of variable 
size is the higher interconnectivity of the grid in the for- 
mer case, which mitigates but not solves completely the 
effects of the biased sampling. The main effect of this 
problem upon our calculations, it is that the size of the 
boundary will be underestimated, and that some parts 
of the landscape are not going to be visited. Thus, the 
method works better before a transition, and in fact only 
provides a bound for the scaling exponents. To improve 
upon this method, some other algorithms can be used, 
as for example collision prediction |38| . 

Now let pk(S) be the probability of state k in phase 
space to be occupied by the system when a grid of scale 
6 is chosen. The random walk process can be viewed as a 
Markov chain j3?J , where the probabilities of visiting each 
microstate are contained in a vector . The probabilities 
at each step are transformed according to an stochastic 
matrix that contains as elements the probability of tran- 
sition among states, 
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where the rows of the matrix are normalized to 1, and 
p k (S) are the probabilities after a step is made. In a 
hard-core system, jumps only occur between allowed grid 
points that are first neighbours. An clement S r t of this 
matrix is zero when one of the states r or t is in the for- 
bidden part of the landscape. S r t is also zero if r or t are 
not first neighbours. The only elements different from 
zero are S r t = l/zt(S), if r and t are allowed neighbours, 
where z t {8) is the coordination in phase space of site t 
(i.e., the number of allowed first neighbours of t) for a 
given scale S. Points at the boundary of the landscape 
are the ones where l/z r (S) < 2DN since they are con- 
nected to points inside the forbidden part of the phase- 
space. A matrix of this type has at least one eigenvector 
with eigenvalue one, while the others have norms equal 
or less than one 9]. Thus, after successive applications 
of the matrix, the stable configuration is given by the 



eigenvalue with norm one,from where it follows that the 
final probabilities satisfy 
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This matrix is similar to the Hamiltonian of a binary al- 
loy in an hypercubic lattice, where the two self-energies 
are very different |39l | (split band regimen). It is easy 
to prove that the final equilibrium vector coincides with 
the bonding state (which corresponds with the maximal 
wave-length of the solution, and nearly zero phase differ- 
ence between sites) of the binary alloy. 

In a Monte-Carlo step, the probability of having a 
rejection is given by the probability of jumping into a 
boundary point (p k (S)), multiplied by the probability of 
jumping form a state k into a state t in the forbidden 
region of the landscape, which is given by the elements 
of the stochastic matrix. Thus, the probability of landing 
in a forbidden state t is, 



( 2DN-z k {8) \ 
PtW={ 2DN ) Pk{6) - 



(6) 



If Cb{5) denotes the set of all boundary points, the total 
probability of having rejections at a scale 8 (we denote 
this probability by p R (5)) is obtained by summing over 
all boundary states k, 



p r (s)= E 

kec B (S) 



Now we write z k (8) as an average plus a fluctuation part, 
Zk(S) =< z k (S) > +Az k (S), where 
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and Mb (6) is simply the number of boundary points at 
a scale 6. The same procedure can be made for Pk(S), 
giving p k (S) = (p B ($)) + Ap k (6), where (pb(S)) is defined 



as. 
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Using these definitions, and that the sum of the fluctua- 
tions is zero, eq.(|7J is rewritten as, 
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The term E] Az k (S)Ap k (S) is a measure of the correla- 
tion between state and coordination fluctuations. Since 
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the eigenvector with eigenvalue one corresponds to a 
bonding state in a binary alloy, using a variational proce- 
dure with a trial function or analyzing the spectral mo- 
ments 39], it can be proved that Apk(S) » Azk(S)/2DN. 
This term gives a correction of order, 
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(10) 



where a (6) is the standard deviation of the phase space 
coordination distribution Zk{8). Thus, it follows that, 
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(11) 



We notice that Mb{5) < pb(8) > is a bound for the 
size of the whole boundary of the landscape. For exam- 
ple, when the sampling is uniform, Mg(8) < Pb(J>) >= 
M B (8)/M(S), since < p B (8) >= l/M{8). The number 

and if the bound- 
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ary has a scaling of the type Mb (5) 
expect a scaling of ea. ljTT|l as, 
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where Df = (D — Db)N is an effective fractal dimension 
due to the different scalings of the grid and the landscape. 
In general, states at the boundary are less visited, thus 
we get the following inequality, 



M B {5)(p B {8)) < 
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(13) 



As a result, we get a bound for this scaling exponent, 
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The evaluation of this bound can be easily implemented 
inside a Monte-Carlo simulation; it only requires the re- 
jection probability p R (8), the average coordination < 
z(8) >, and the fluctuations <5(<5). To do this, first we 
divide the phase space with a grid of spacing S. Then 
we perform the Monte-Carlo simulation, but if there is a 
rejection during a trial step, this means that the initial 
state is in the boundary of the landscape. To look at the 
coordination in phase space of this state, a movement in 
each of the DN coordinates of the gird is performed, as 
in eq.©, but for each coordinate q r taken in sequence 
from r — 1 to r = DN. For each coordinate move- 
ment, the new state is tested in order to determine if 
its an allowed or forbidden state. After this cycle in the 
coordinates, the number of accepted states is the coordi- 
nation number Zfe(<5). The process continues until a new 
rejection appears, and at the end of the simulation, the 



FIG. 1: Two disks in a rectangular box of length L and width 
a. Below the box, the corresponding configurational part of 
the phase space is shown. The allowed part of the landscape 
is the area indicated with the grey shadow. A grid of scale 
8 is indicated by dotted lines. Boundary points are indicated 
by open circles. Closed circles are states in the forbidden part 
of the phase-space. Notice that in this problem, ergodicity is 
always broken, since the allowed parts of the landscape are 
not connected. 



average and the standard deviation of the distribution of 
Zk(8) are obtained. The same procedure is repeated for 
different scales S. 

Figure 1 illustrates the procedure for a very simple sys- 
tem. Consider two disks in a box of length L and width 
(j. In such a case, the movement is one dimensional. 
If x\ and 22 are the coordinates of each disk, the con- 
figurational part of the phase space can be represented 
as a plane. The shape of the landscape is determined 
by the condition of non-overlap \x\ — Xi\ > a and the 
walls of the box. The allowed phase space is made by 
two triangles, as shown in fig. 1. Notice that ergodic- 
ity is always broken since the two triangular regions are 
not connected. In fig. 1, the grid is indicated as dotted 
lines; the points at the boundary (open circles in fig. 1) 
are those connected to grid points that are outside the 
triangular regions (closed circles). 

In this simple system, it is very instructive to compare 
the rejection ratio of the Monte-Carlo with the topol- 
ogy of the phase space, since the theoretical value for 
M B (8)/M(8) for a given packing fraction <fr = 2/L is 
easy to find. The value of Mb{8)/M{8) in this case is 
given by the ratio between perimeter-area of the triangle 
as a function of the scale, 



M B {6) 
M{8) 



2(2 + y/2) 
L(l-2a/L) 



where the subscript is used to indicate that the result 
depends on the corresponding length of the box. This 
result can be related with the probability of rejection of 
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FIG. 2: Rejection ratio as a function of the scale <5 (mea- 
sured in units of a) for L = 50a (squares) , 30<r (circles), 15a 
(diamonds), 10a (triangles up) and 4a (triangles down). 



FIG. 3: Slopes of the fitting lines that appear in fig. 2 as a 
function of L. The solid line is the function m(L). 



the Monte-Carlo as follows. If we suppose an uniform 
sampling, the probability of hitting a boundary point is 
given by the perimeter- area ratio of the triangles. The 
average coordination of the boundary points can be ob- 
tained from a direct inspection of the drawings with dif- 
ferent grids, that gives (z(S)) ~ \/2 + (3/2) for a < L. 
We can neglect the term <t(o)/4, since it is very small to 
be considered (this approximation was confirmed after- 
wards by the coordination statistics obtained from the 
Monte-Carlo simulation). Using eq. (|llf) . the predicted 
rejection probability is, 

p R (S) = m(L)S, 



m(L) is well approximated by, 

m (L)«*£±^(l + 2a/L), 

so m(L) r~j L^ 1 , and Df = 1 as expected for a smooth 
surface. In the region ctrjL, the difference between both 
results is due to the fact that the average coordination 
number is not anymore \[2 + (3/2), and a correction is 
needed in the analytical formula. Also, in this region the 
sampling is far from uniform, since the grid is very small 
compared with the size of the boundary. 



where m{L) is defined as, 

m(L) = ^-^=)[M B (6)/M(S)} L . 

The rejection is thus expected to be proportional to S, as 
confirmed in Figure 2 by the numerical simulations using 
a Monte-Carlo, where the rejections are plotted against 
5 for different L. 

Using a least-square fitting, the slopes for each of the 
lines is shown in figure 3 using a log-log plot. The solid 
curve is the theoretical value of m(L) and the squares 
are the results of the simulation using the Monte-Carlo 
simulation. These results are in very good agreement 
with the theoretical values, specially for a <C L, where 



III. SCALING IN A SIMPLE LIQUID 

In this section, we show the results obtained using 
the method proposed in the previous section for a two- 
dimensional system composed of hard-disks at different 
packing fractions <j> = Nira 2 /AS, where N = 100 parti- 
cles. The thermodynamics of this system has been in- 
vestigated by many different groups during the last 50 
years [4(j [4l]] [43 ■ Despite the simplicity of the model, 
the nature of the phase transition from solid to fluid 
is still debated |43|. as well as the nature of local or- 
der 0] and its relation with some peaks in the radial 
distribution function |45j . Here we will only investigate 
the landscape scaling. Figure 4 shows a log-log plot of 
Mb{8) < pb(S) > as a function of 5 for different packing 
fractions, as indicated in the figure caption. 
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FIG. 4: Parameter Mb (S) {pb(S)} as a function of the scale S, 
for different packing fractions. From top to bottom, <f> = 0.74 
(squares), </> = 0.71 (x),(f> = 0.39 (triangles), <f> — 0.23 (stars), 
tf> = 0.12 (diamonds), 4> = 0.08 (squares) and <j> = 0.04 (filled 
circles). The lines were obtained using a power law fit. 



Figure 5 is a similar plot, but only for packing frac- 
tions near the freezing point (denoted by <p ). Both 
plots give evidence that there is a power law scaling of 
Mb (5) < pb(S) > with S. This power law behavior is 
clear near the freezing point or at low densities, where 
fits of the type S Df are shown for the different sets of 
data (only fits with correlation coefficients bigger than 
99% are shown). Notice that all the fits have a cut-off 
at S = 0.05<7, since there is a cross-over in the power law 
behavior, i.e., for a given packing fraction, two regions 
with different scaling exponents are observed, as seen in 
Figure 5, where a drop of Mb{$) < Pb{$) > is observed 
around 5 = 0.05o\ For low packing fractions, the expo- 
nents for 5 < 0.05cr tend to be smaller than in the region 
5 > 0.05<7. The fact that two exponents are observed 
means that below a certain length-scale, the landscape 
has a different structure. For all the different packing 
fractions, this behavior is nearly similar. We can spec- 
ulate that this change of regimen for the scaling at a 
length-scale is related with the different processes of re- 
laxation that have been observed in diverse simulations 
[28| |4^ | fl?! and experiments ^(| , since although a Monte- 
Carlo simulation does not provide the real dynamic of the 
system, is clear that a big length scale 6 in phase space 
correspond to long times in the evolution of the system, 
as also expected from the Adam-Gibbs relation between 
relaxation times and configurational entropy 0. How- 
ever, this speculation needs to be investigated in more 
detail. 

We also notice that for packing fractions 0.2 < (f>/(f> < 
0.6, it seems that using one single scaling exponent is 
not enough to fit the data, which is an indicative of a 
multifractal structure, although if we restrict the fitting 
for 6 > 0.2, again a good power law fit is obtained. 



FIG. 5: Parameter Mb{5) (pb(S)) close to freezing as a func- 
tion of S, at <f> — 0.74 (squares) and <f> = 0.63 (circles). 



In figure 6 we plot the scaling exponents obtained from 
the data of fig. 4 and 5 as a function of the packing frac- 
tion, for the regions (0.2cr < 5) where a clear scaling is 
obtained for all the graphs . As shown in the figure, 
as the packing fraction reaches the freezing point, Df 
goes to zero, and the landscape boundary scales nearly 
as the volume in phase space. This means that near the 
freezing point, the topology of the landscape restricts in 
a severe manner the available phase space. Thus, fig- 
ure 6 provides clear evidence of how the topology of the 
landscape is responsible for the phase transition that oc- 
curs at the freezing point, and reinforces the speculation 
about relaxation times, since it has been observed in ex- 
periments with colloids that freezing occurs when long 
time relaxation is not longer available [l6| . 



IV. CONCLUSIONS 

In this article, we have discussed some aspects of how 
to characterize the structure and texture of the energy 
landscape in simple fluids. As a result, we showed a 
method to investigate the boundary of the landscape that 
uses the Monte-Carlo rejection ratio plus the average co- 
ordination of a state in phase space. An example of how 
to apply the method has been presented for a very sim- 
ple model that consists of two disks that moves in one 
dimension. A similar procedure applied to a system of 
hard-disks shows a clear power law scaling of the ratio 
between the boundary and the volume of the landscape. 
A cross-over in the scaling exponents has been observed 
for a given packing fraction. Near the freezing point, the 
boundary of the landscape scales as the volume in phase 
space, and as a result the system tend to stay in pock- 
ets of the phase-space. We speculate that the cross-over 
observed in the scaling is related with the different kinds 
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of relaxation processes of the fluid. In future works, we 
will further explore this idea. 



FIG. 6: Exponent Df as a function of the ratio between the 
packing fraction and the packing fraction at freezing (<j> ). 
The size of the squares is proportional to the maximal error, 
and the line is a visual guide to the eye. 
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